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Abstract 

The convergence of the iterative solutions of the transport equations of cosmic muon 
and tau neutrinos propagating through Earth is studied and analyzed. For achieving 
a fast convergence of the iterative solutions of the coupled transport equations of v T , 
v T and the associated fluxes, a new semi-analytic input algorithm is presented 
where the peculiar r-decay contributions are implemented already in the initial 
zeroth order input. Furthermore, the common single transport equation for muon 
neutrinos is generalized by taking into account the contributions of secondary v„ 
and Ufj, fluxes due to the prompt r-decay r — ► initiated by the associated tau flux. 
Differential and total nadir angle integrated upward-going \i~ + fi + event rates are 
presented for underground neutrino telescopes and compared with the muon rates 
initiated by the primary z/,j, v t and r fluxes. 



1 Introduction 



Upward-going cosmic neutrinos with energies below 10 8 GeV play a decisive role for un- 
derground neutrino telescopes, since the atmospheric background can be more effectively 
controlled, in contrast to downward-going cosmic neutrinos. While traversing through 
the Earth, upward-going muon (anti) neutrinos undergo attenuation (absorption) due to 
weak charged current (CC) and neutral current (NC) interactions as well as regeneration 
[f , 2] due to NC interactions. The latter shift the energy of the neutrinos, rather than 
absorbing them, to lower energies and populate the lower energy part of the initial cosmic 
neutrino flux spectra, thus adding to the naive non-regenerated fj,~ + fj, + event rates at 
the detector. Such propagation effects of muon (and electron) neutrinos through Earth 
are described by a single transport (integro-differential) equation which can be rather 
easily solved iteratively [1, 2, 3, 4, 5]. 

On the other hand tau (anti)neutrinos are not absorbed, but degraded in energy, in 
the Earth as long as the interaction length of the produced tau leptons is larger than 
their decay length (which holds for energies up to about 10 9 GeV). Because of these latter 
(semi)leptonic decays r — > v T X, the Earth will not become opaque to v T [6] since the r~ 
produced in CC interactions decays back to v T . This 'regeneration chain' v T — > r — > u T — > 
. . . continues until the v T and v T , as well as the r ± leptons, reach the detector on the 
opposite side of the Earth. Thus the propagation of high-energy tau neutrinos through 
the Earth is very different from muon and electron neutrinos, and we have now to deal 
with coupled transport equations for the T and fluxes [4, 7, 8, 9, 10, 11, 12, 13]. 
Obtaining stable iterative solutions of these coupled integro-differential equations is far 
more involved as compared to the single transport equation for muon neutrinos. It is 
one of our main objectives to discuss the general qualitative and quantitative structure 
of these solutions and to present an efficient input algorithm which allows for a rather 
fast convergence of the iterative procedure. This applies to all present model cosmic 
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neutrino fluxes. Moreover the r~ + t + flux, generated by the initial cosmic v T + v T flux 
while traversing the Earth, gives rise to a secondary + flux [14] via r — > z/ M due 
to the prompt r-decays like r~ — > v T pTv >1 . This adds considerable contributions to the 
primary cosmic + flux and may increase the fj,~ + fi + rates at the detector site 
sizeably [9, 12], depending on the cosmic flux and nadir angle considered. Such effects 
require an extension of the simple single transport equation for and the inclusion of 
the appropriate prompt decay term reduces the convergence of the iterative procedure 
considerably. 

The simple single transport equation for ^ will be discussed for completeness in 
Sec. 2. Although frequently used, the excellent convergence of its iterative solutions 
has not been explicitly demonstrated thus far for more realistic cosmic neutrino fluxes, 
apart from some specific steep toy model neutrino fluxes [2]. In Sec. 3 we turn to the 
iterative solutions of the far more complicated coupled transport equations for T and 
their associated fluxes. A new semi-analytic input algorithm is presented which allows 
for a fast convergence of the iterative solutions. The implications for the upward-going 
pT + /i + event rates for underground neutrino detectors for some relevant cosmic neutrino 
fluxes will be briefly outlined as well. The solutions of the generalized single transport 
equation for muon neutrinos, by taking into account the contributions of the secondary 
v n + fl ux from prompt r 1 * 1 decays based on the calculated associated r 1 * 1 fluxes, are 
discussed in Sec. 4. Their implications for the expected \i~ + /i + event rates, as initiated 
by various relevant cosmic neutrino model fluxes, are presented as well. Finally, our 
conclusions are summarized in Sec. 5. 



2 The transport equation of muon neutrinos 

Disregarding possible contributions from other neutrino flavors for the time being, the 
transport equation for upward-going cosmic muon (anti) neutrinos ^ passing through 
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Earth can be written as [1, 2, 3, 4, 5] 

dF u (E,X) F V (E,X) , 1 f 1 dy 



f J^K™(E,y)F^E y ,X) (1) 

Jo 1 y 



OX X U (E) \v(E)J 1 

where F v = d<& u ^/dE is the differential cosmic neutrino flux and E y = E/(l — y). The 
column depth X = X(9), being the thickness of matter traversed by the upgoing leptons, 
depends on the nadir angle of the incident neutrino beam {9 = 0° corresponds to a beam 
traversing the diameter of the Earth); it is obtained from integrating the density p(r) 
of the Earth along the neutrino beam path L' at a given 9, X(9) = p(L')dL' with 
L = 2_R e cos 9, R e ~ 6371 km, denoting the position of the underground detector, and 
X(9) can be found, for example, in Fig. 15 of [15] in units of g/cm 2 = cm we. Furthermore 
A" 1 = N A al%, Na = 6.022 x 10 23 g -1 , is the inverse neutrino interaction length where 



_tot _ CC , NC J 



1 da™(E s , V ) 

The various CC and NC ^uN cross sections are calculated as in [5, 13], with the relevant 
details to be found in [16], utilizing the QCD inspired dynamical small-x predictions for 
parton distributions according to the radiative parton model [17]. Notice that conven- 
tionally fitted parton distributions at the relevant weak scale Q 2 = would require 
additional ad hoc assumptions (see, e.g., [15, 18]) for the necessary extrapolations into 
the yet unmeasured small Bjorken-x region x < 1CT 3 (x ~ M^/2m N E). The first term 
in (1) describes the attenuation (absorption) of neutrinos when penetrating through the 
Earth, and the second one their regeneration consisting of the degrading shift in their 
energy. For definiteness all formulae are given for an incoming neutrino beam, but similar 
expressions hold of course for antineutrinos. 



Equation (1) can be efficiently solved by the ansatz [2] 

X 



F^(E,X)=F^(E)exp 



(3) 
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with an effective absorption (interaction) length 

and where F^(E) = F Vli (E, X = 0) denotes the initial cosmic neutrino flux which reaches 
the Earth's surface. Depending on the assumed cosmic neutrino flux, the Z^-factor can 
take any non-negative values. Its physics interpretation and the consequences for the 
shadowing factor S = exp [—X/A u ] in (3) are immediate: Z v < 1 (the only case considered 
in [2] relevant for steeper, i.e., soft model fluxes) implies A„ > \ v > thus S < 1, 
i.e. the neutrino flux will be further attenuated since absorption plays the dominant 
role; for Z u — 1, A u — oo, i.e. S = 1 which means that regeneration and absorption 
compensate each other; finally Z u > 1 implies A u < and S > 1, and consequently the 
NC regeneration in (1) can even cause an enhancement of the neutrino spectrum with 
respect to the initial flux F^(E) for certain energies and depths X. Inserting (3) into (1) 
yields 

Z^(E,X) = 1 jf'iY' jf^y^ c (S,y)^(^y) e -^(W') (5) 

with Vu (E,y) = F^(E y )/(l - y)Fl(E) and D^(E,E y ,X>) = A^(E y ,X') - A^(E,X'). 
Using an iteration algorithm to solve for Z U/1 (E, X), one can formally rewrite the solution 
of (5) after the n-th iteration as 



(6) 



where 

n (n)/p p va 1 - Z Ufi (E y ,X') 1 - Z Ufi (E,X') 

D W E > E v> X ) = TJe) • (7) 

The reason why this iteration is expected to converge very fast is as follows: the kernel 

K^ c peaks very strongly [2, 19] at y = and y — 1, with the contribution at y ~ 1 

being, however, exponentially suppressed in (6); thus the main contribution to the integral 

over y in (6) comes from the region around y ~ where D U ^(E, E y , X') — > as y — > 0. 

Therefore the iteration should be robust with respect to choosing the n = approximation 
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[2]. The most simple input choice is ZlJ(E, X') = in (7). For this case the analytic 
X'-integration in (6) yields 

f 1 1 _ p -XD v {E,E y ) 

Z^(E,X) = j( dyK^°(E, y) r) v (E, y) YdJe^) (8) 

with 

DJE.EJ^D^E^X^^-^. (9) 

With the n — 1 solution in (8) at hand, it is now straightforward to obtain iterations in 
higher orders, for example, for n = 2 by inserting (8) into (7) gives zff in (6). 

Representative cosmic neutrino fluxes of some hypothesized sources are displayed in 
Fig. 1 which we shall partly use for all our subsequent calculations. Recent diffuse neutrino 
flux upper limits of AMANDA [20, 21] are shown by the bars with arrows - the latter 
indicate the still allowed region. Although the huge flux from active galactic nuclei of 
Stecker and Salamon (AGN-SS) [22] has been already excluded, we shall use it merely 
as a theoretical playground due to its unique spectrum at lower energies where (E) ~ 
const, for E < 10 5 GeV. On the other hand the AGN-M95 flux [23] is still compatible 
(although slightly in conflict) with the AMANDA upper bound, as are the gamma ray 
burst (GRB-WB) [24] and topological defect (TD-SLBY) [25] fluxes. These latter three 
fluxes will be used for our 'realistic' model calculations. The TD-SLSC [26] and Z-burst 
[27] fluxes are shown just for illustration since they are too minute for being tested with 
upward-going event rates [5]. Note that the initial cosmic (anti) neutrino fluxes F® P (E) 
in (3) which reach the Earth's surface are given by F^ = = F° t = F? t = \d$/dE 
with $ being the cosmic + flux at the production site in Fig. 1. 

For a better comparison of our quantitative results with the ones obtained in the 
literature, we also employ two generic initial fluxes incident on the surface of the Earth 
at a nadir angle 6 = 0° of the form [4, 7] 

F^(E) = N.E-'il + E/Eo)- 2 , £ = 10 8 GeV (10) 

= N * E ~ 2 (ii) 
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with adjustable normalization factors N, h for example, N± — \ x 10 13 /(cm 2 sr s) and 
N 2 = \ x 10~ 7 GeV/(cm 2 sr s). Notice that the generic E 1 energy dependence is 
representative for the TD and Z-burst fluxes in Fig. 1 for E 10 7 GeV; and also for the 
GRB-WB flux for E < 10 5 GeV. Furthermore the latter GRB-WB flux behaves like E' 2 
in (11) for 10 5 < E 10 7 GeV, where such a power spectrum with index —2 is typical 
for shock acceleration (see, e.g., [21]). 

Our results for Z$ and zf?J are shown in Figs. 2 and 3 for two typical values of the 
nadir angle, 9 = 0° (X = 1.1 x 10 10 cm we) and 9 = 50° (X = 3.6 x 10 9 cm we). The 
iteration converges very fast since in general the maximum difference between Z^J and zj?J 
is less than about 5%, \zj?J /Z^J — 1| < 0.05, and moreover \Z^J /Z$?J — 1| < 0.005. Thus 
the first n — 1 iteration is already sufficiently stable and suffices for all cosmic neutrino 
fluxes considered at present [19]. Notice that for larger 9 (smaller X) the difference 
between zj?J and Z^J decreases and therefore the stability increases. The results for Z p ^ 
are similar but Z P/i > Z Vfi for E < 10 6 GeV where \ p > \ u . The resulting and fluxes 
follow from (3) and can be found in [4, 5, 7]. 



3 The transport equations of tau neutrinos and taus 



Apart from the absorption (attenuation) due to <t*% and regeneration due to 0^$ in (1), 
for upward-going cosmic tau (anti) neutrinos ^v T , it is important to take into account the 
regeneration from the t^ 1 decays as well as the contributions from the CC tau interactions. 
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The tau neutrino and tau fluxes then satisfy the following coupled transport equations: 
dF Vr {E,X) F Ut (E,X) 1 f 1 dy NC 

+ / -^-K T (E,y)F T (E y ,X) (12) 

Jo 1 — y 

dF T (E,X) F T (E,X) d[7(E)F T (E,X)] 

dX \(E) dE 

+ X^l!^y K ^ y)F ^ X) (13) 

where F Vt = d& Ur /dE and F T = d& T /dE are the differential energy spectra (fluxes) of tau 
(anti) neutrinos and tau leptons and the initial fluxes at the surface of the Earth (X = 0) 
being given by F® t (E) = F^(E) = \d<b/dE with $ being the + flux at the cosmic 
production site in Fig. 1. The cross section kernel K^ c is defined in (2) and a similar 
expression holds for K^ c . Furthermore 

K r(E, y) = ^ K? C (E, y) + K^(E, y) (14) 

where 

K CC (F v) - 1 d(J rS(E y ,y) d _ 1 dT T ^ UTX ,(E y ,y) 

and A" 1 = N A a™ = N A (a™ + a™), and A" 1 = (A^ c ) _1 + (A?*)" 1 with (A^ c ) _1 = 
N A aCC in (13). The decay length of the r± is \* CC (E,X,9) = (E/m T )cT T p(X,9) with 
m T = 1.777 GeV, cr T = 87.11 jum and p denoting the Earth's density (see, e.g., [15]). 
Furthermore, since l/T t ° t (E) = (E/m T )r T , the r-decay distribution in (14) becomes 
K* ec {E, y) = (l- y) dn{z)/dy with z = E U JE T = E / E y = 1 - y and [7, 28] 

^ = E s *boW+^iw] as) 

with the polarization P = ±1 of the decaying r ± . The r — > z/ r X' branching fractions 
_Bj into the decay channel i and the functions g^^z) are given in Table I of [7]. The 



decay channels i considered are r — > v T [iv^ r — > z/ T 7r, r — > z/ T p, r — > z/ r ai and r — > z/ T X 
which have branching fractions of 0.18, 0.11, 0.26, 0.13 and 0.13, respectively. The lepton 
energy-loss is treated continuously [29, 30, 31] by the term proportional to 'y(E) in (13). 
Alternatively, the average energy-loss can be treated separately (stochastically) [32, 33], 
i.e., not including the term proportional to j(E) in (13) but using instead —dE/dX = 
'y(E) = a + f3E. We shall compare these two approaches for taus and muons toward the 
end of this Section. The most general solution of Eqs. (12) and (13) has been presented 
in [10, 13], and in the context of atmospheric muons in [31]. For the time being, however, 
we disregard the 7-term in (13) since observable non-negligible upward-going event rates 
are obtained only for energies E < 10 8 GeV [7, 13] where the energy-loss of the taus can 
be neglected [10, 32, 33, 34, 35]. 

In the relevant energy region below 10 8 GeV, the tau-lepton interaction length is much 
larger than the decay length of the r (see, e.g., [33] and below), X T (E) 3> Xf ec (E), i.e., 
K T ~ K^/Xf™ in (14) and A; 1 ~ (X^)- 1 in (13). Solving (12) and (13) with a similar 
ansatz as for muon neutrinos in (3), we write 

X 



F l/T (E,X)=F° T (E)exp 
with an effective interaction (absorption) length 



(16) 



^ x ^Y^x) ^ 

where Z = Z Vt + Z T . Inserting (16) into (12) and (13) yields [4, 13] 

Z Ut (E,X) = ^ £ dX' j\yK^ c (E,y)r] v (E,y)e- x,D ^ E y' x ') (18) 
with Vu as in (5) since F^ = F° and D Vr (E, E y , X') = A^(E y , X') - A^(E, X'), and 
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where the obvious dependence of A^! cc on 9' = O(X') has been suppressed and 




'yy'' 



(20) 



X" 



with Eyyi = Ey/ (1 — y') — E / ( 1 — y) (1 — %)'). Notice that the tau-flux F T is generated by 
the CC interactions of the initial F° t flux and attenuated in addition due to its decay. In 
order to solve for Z(E, X) iteratively as for the fluxes in the previous Section, one has 
to make a proper choice for the initial input. Due to the D Vt function in the exponential in 
(18) with D Vt — > in the relevant y — * region, the iterative result for Z Vt {E, X) is very 
robust with respect to the initial input choice, as discussed after (7). Therefore we use 
again Z$(E, X') = on the rhs of (18). In the case of Z T (E, X) in (19) there is, however, 
no equivalent exponential as in (18) and thus the convergence of the iterative procedure 
becomes sensitive to the input choice. It turns out that a convenient and efficient input 
choice can be obtained by implementing the peculiar E and X dependence as implied 
by the r-decay contributions in (19) from the very beginning. This can be achieved by 
choosing a vanishing Z-factor on the rhs of Z T in (19), in which case the X'-integral can 
be performed analytically [36] and the input for the total Z-factor becomes [13] 



where the last approximation is due to A!; cc <^ \ v in the relevant energy region E < 10 
GeV, i.e., Z^ becomes practically independent of the decay length A^! ec . Furthermore 



Z<®(E,X) 




J 1 (l _ e -XD TV {E,E y ) 

XD VT [E yi Eyyi) lD T „(E,E y ) 
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D U (E, E y ) is given in (9), D UT (E,E y ) = l/\ v {E y ) - l/\f c (E,6) and D TU (E,E y ) = 
—D L , T (E y , E). We have checked that this input guarantees, for all cosmic neutrino fluxes 
considered at present, a faster convergence of the iterations than choosing [4] the solution 
for the flux as an input, = Zg } with Z$ given in (8). Moreover, choosing [10] 
a vanishing initial input, Z^> = 0, as was perfectly sufficient for the fluxes, results in 
the worst, i.e., slowest convergence of the iterative procedure. One can now rewrite the 
solution for Z(E,X) in (18) and (19) after the n-th iteration as [36] 



Z^(E,X) = 1 J* dX' j\yK^(E,y) Vv (E,y)e- x ' D ^ E ^ 

+ x^(e,o) Jc £ dx ' I! dy K * ec{E ' y) vAE > y ^ Kl ^ 



xe -X'/X^(E y ,e) e X'/A^(E,X') f X dx n [' dy > K CC yl) y/) 

Jo Jo 

x e -X"/A^{E yyl ,X") e X"/4<*{E y ,8) ( 2 2) 



where A£\E,X') = X V (E)/[1 - Z^(E,X% i.e., 

D ^(E,E y: X)- ^ ^ . (23) 

Accordingly, the iterations have to be started with our initial n = input in (21). After 
having obtained the final convergent result for Z^ n+1 \ the final v T flux F^ +1 \E, X) follows 
from (16), 

F^ +1 \E,X) = i^ (jB)e -x/A£ +1) (£,x) , (24) 



which in turn gives the r-flux 

F^ n+1 \E, X) = — l — e "WW) f X dX' f K° C (E, y) F<£ +1 \E y , X') e x '/ x ^ E ^ 
K{E) J J 1 — y 

(25) 

Similar expressions hold for antineutrinos as well. 

The iterative results for the total Z-factor in (17) are shown in Figs. 4 and 5 where 
the initial input Z^°\ as given in (21), is displayed by the dotted curves. For the generic 
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initial E~ l and E~ 2 fluxes in (10) and (11) we also show in Fig. 4 the results after the 
third iteration, Z^\ in order to illustrate the rate of convergence as well as its dependence 
on the nadir angle 6 = 0° (X = 1.1 x 10 5 km we) and 6 = 50° (X = 3.6 x 10 4 km we). In 
general it turns out that already the second (n = 2) iteration yields sufficiently accurate 
results, Z^ 2 \ provided one uses as input Z^ in (21) as implied by the r-decay. This holds 
also for the rather hard initial E~ x flux in Fig. 4 and the AGN-SS flux in Fig. 5 which 
imply large Z-factors, Z>1. This is so because the maximum difference between the 
results of the next n = 3 iteration Z^> and Z^ is less than about 5% for all relevant initial 
cosmic neutrino fluxes. An accuracy of less than about 5% is certainly sufficient in view 
of the uncertainties inherent to models of cosmic neutrino fluxes (cf. Fig. 1). Obviously 
the iterative convergence improves even more for larger values of 9, i.e., smaller depths 
X, as can be deduced from Fig. 4. It should be emphasized that, in contrast to the case 
of muon neutrinos in Sec. 2, the first n = 1 iterative results for Z^ are not sufficiently 
accurate as can be seen from Figs. 4 and 5 by comparing the dashed curves (Z^) with 
the solid ones (Z^): in some cases (harder initial fluxes) Z^ becomes larger than Z^ 
by about 20%. The results for P T , Z = Z Pt + Z T +, are again similar but larger than for u T , 
Z = Z Vr + Z T - , for E <, 10 6 GeV where \ p > A„. Inserting the various iterative solutions 
of Figs. 4 and 5 into (16) we obtain the v T fluxes for a given n-th iteration, F^\E,X). 
The ratios of these fluxes for two consecutive iterations, F^ +V> / are displayed in 
Figs. 6 and 7. Whereas the first iteration relative to the zeroth input, F$ /F$, is way 
off the final result as shown by the dashed curves, the second iteration suffices already for 
obtaining a sufficiently accurate result as illustrated by F$ ' / F$ by the solid curves. This 
is supported by the fact that an additional third iteration results in \F$ /F$ — 1| < 0.05 
for all relevant initial cosmic neutrino fluxes considered. (This stability does not hold [7] 
for initial fluxes F® +D ~ E^ 1 without an appropriate E~ 2 cutoff in (10) at very high 
energies, or for fluxes which are partly even flatter than E^ 1 up to highest energies like 
the Z-burst flux in Fig. 1. This instability is caused by the fact that i] u (E,y) = 1 for 
F° ~ E' 1 in Z T in (19); thus the huge spike of da^/dy at y -> 1 in (19) and (20) does 



11 



not get damped by powers of (1 — y) - as opposed to, e.g., F° ~ E~ 2 where r\ v — 1 — y. 
This, however, is of no concern for Z v and Z Vt in (5) and (18), respectively, since there 
the integrands are exponentially suppressed as y — > 1 via exp [— X'D V (E, E y , X)]. Notice 
that the Z-burst flux is far too small for being tested with upward-going muon events 
[5, 13].) Therefore we consider F ( (2 J (E, X) as our final result. It is furthermore obvious 
from Figs. 6 and 7 that the convergence of the iterative procedure strongly improves for 
increasing values of 9 (decreasing X) as illustrated for 9 = 50°. 

The resulting total v T +P T fluxes for the various initial total cosmic fluxes F® t+Pt (E) are 
shown in Figs. 8 and 9 for three typical nadir angles 9 — 0° (X — 1.1 x 10 5 km we), 9 = 30° 
(X = 6.8 x 10 4 km we) and 9 = 60° (X = 2.6 x 10 4 km we). The typical enhancement 
('bump') of the attenuated and regenerated ^v T flux around 10 4 - 10 5 GeV at small values 
of 9, which is prominent for harder (flatter) initial fluxes like ~ E~ l in Fig. 8, and 

which is absent for v M fluxes, agrees with the original results of [4, 7, 9] , as was also 
confirmed by a Monte Carlo simulation [34]. Such an enhancement is less pronounced 
for the GRB-WB and TD-SLBY fluxes in Fig. 9, and is absent for steeper fluxes like for 
the E~ 2 one in Fig. 8 and for the even steeper AGN-M95 flux in Fig. 9. Regeneration is 
responsible for an even more pronounced enhancement below 10 4 GeV for the AGN-SS 
flux in Fig. 9 since this flux is particularly hard below 10 5 GeV (cf. Fig. 1). This latter 
result is shown mainly for theoretical curiosity. From now on we shall disregard the cosmic 
AGN-SS flux since it is in serious conflict with recent experimental upper bounds [20, 21] 
as can be seen in Fig. 1. The results for the absolute total v T + v T and t~ + r + fluxes, 
arising from the initial cosmic v T + u T fluxes, are presented in Figs. 10 and 11. The u T + u T 
results correspond of course to the relative ratios shown in Figs. 8 and 9. Besides the 
generic initial fluxes in (10) and (11), we have in addition used only those initial cosmic 
fluxes in Fig. 1 which give rise to large enough upward-going fi~ + fi + event rates [5, 13] 
measurable in present and future experiments. Note that the r~ + r + fluxes in Figs. 10 
and 11 at the detector site, despite being (superficially) suppressed with respect to the 
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v T + v T fluxes, sizeably contribute to the upward-going fi~ + /i + and shower event rates 
[13]. This is due to the fact that the r fluxes do not require additional weak interactions 
for producing /i-events in contrast to the v T fluxes. Because of the prompt decays, 
they furthermore give rise to a sizeable secondary + flux contribution to the original 
cosmic + Dp flux, which will be discussed in Sec. 4. 

In [13] a semi-analytic solution of the coupled transport equations (12) and (13) has 
been presented and used, which was obtained from the first n — 1 iteration starting with 
a vanishing input Z°(E,X) = 0, instead of using (21). As we have seen, this approach 
does not provide sufficiently accurate results, despite opposite claims in the literature [10] 
(the first n — 1 iteration is sufficient only for very large values of 9 close to 90°, i.e., very 
small values of X/p = (9(100 km), relevant for neutrinos skimming the Earth's crust). 
This n — 1 iterative solution of [13] underestimates the correct results in some extreme 
cases (like for the hard initial E^ 1 flux at 9 — 0°) by as much as 40%. On the other 
hand, for increasing values of 9 this discrepancy disappears very quickly. Consequently, 
some of the total nadir-angle-integrated upward-going \i~ + p + event rates calculated 
in [13] will be increased by less than about 2%. This is due to the fact that 80% of the 
pT + p + rates are initiated by the + flux and only about 20% derives from the u T + V T 
and the associated t~ + r+ fluxes. For completeness we present in Table 1 the correct 
expectations for the total fx~ + fi + event rates for the relevant dominant initial cosmic 
fluxes in Fig. 1, using Eqs. (12) and (14) of [13] for calculating the rates initiated by the 
v T + v T and r~ + r + fluxes, respectively. 

Finally, it is also of interest to compare the tau-lepton range as given by our semi- 
analytic approach of treating the energy loss continuously in (13), with the one obtained 
by a stochastic treatment of the lepton energy loss (where the 'j(E) term in (13) is absent, 
i.e., the energy loss is treated separately, and the relevant survival probability P(E,X) 
is calculated using Monte Carlo simulations, e.g., [32, 35, 37, 38]). To do this, we can 
drop the inhomogeneous neutrino term in (13) and the resulting homogeneous transport 
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equation for F T (E,X) can be easily solved [13]: 
F T (E, X) = F T (E(X, E), 0) exp 



x 



A{E{X',E))dX' 



(26) 



with A{E) = 1/\{E) - fry(E)/dE, A" 1 = (A^)- 1 + (A^ )" 1 , and where dE(X, E)/dX = 
^(E) with E(0, E) = E. The survival probability P(E , X) for a tau-lepton with an initial 
energy E at X = is then defined by the ratio of the energy-integrated differential 
fluxes F T at X and X = : assuming, as usual, a monoenergetic initial flux in (26), 
F T {E{X, E), 0) ~ <J(E - E ), one obtains [10] 



p(£ °' A ' )= ^t! exp 



A(£b(X , , J E*,))dX' 



(27) 



where we have used [13] dE/dE = j(E)/^(E) and dE (X,E )/dX = -j(E ) with 
Eq(0,Eq) = Eq. The (tau) lepton range for an incident lepton energy E and a final 
energy E(X, E) required to be greater than E mm at the detector, say, is then defined by 

R(E) = / P(E,X)dX (28) 
Jo 

where we have substituted E for Eq in (27) and the upper limit of integration X max 
derives from E(X, E) > E mm . (Notice that for energy-independent values of a and (3 
in 'y(E) = a + (3E one simply gets X max = -| In a "^f in .) For calculating the r-lepton 
range R T (E) we use in 7 T for the ionization energy loss [39, 40] a T ~ 2.0 x 10~ 3 GeV 
(cm we)" 1 and for the radiative energy loss (through bremsstrahlung, pair production 
and photonuclear interactions) [10] f3 T = f3 T {E) ~ [0.16 + 0.6(£/10 9 GeV) - 2 ] x 10~ 6 (cm 
we) -1 which parametrizes explicit model calculations [32, 35] for standard rock (p = 2.65 
g/ cm 3 ) reasonably well for 10 3 £ E £ 10 9 GeV. Furthermore we impose [32] E mm = 50 
GeV. Our results for the r-lepton range are shown in Fig. 12 which agree of course with 
the ones in [10]. The r-decay term dictates the r-range until E > 10 7 GeV where the 
tau-lepton energy loss becomes relevant. The dashed-dotted curve shows the range as 
obtained by omitting the contribution due to the CC interaction length A^? c in A in (27). 
This term will be relevant for E £ 10 10 GeV where A^ c becomes comparable to A^ cc as 
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evident from Fig. 13. For comparison, stochastic Monte Carlo evaluations [32, 35, 38] 
of the r-range are shown in Fig. 12 by the dotted curves which are of course strongly 
dependent on the assumed model extrapolations of the radiative cross sections to ultrahigh 
energies. Our results depend obviously also on such extrapolations due to specific choices 
of f3 T (E). Nevertheless, one concludes [10] from Fig. 12 that the continuous tau-lepton 
energy loss approach, as used in (13), yields very similar results as the stochastic Monte 
Carlo calculations where the energy loss is treated separately. 

A similar conclusion holds for the muon-range R^(E) which we show for completeness 
in Fig. 14. Within the continuous muon energy loss approach, follows from (28) 
and (27) where the A^ cc term has to be omitted and in 7 (U (-E') = + fl^E we take 
a M ~ a T ~ 2.0 x 10~ 3 GeV (cm we)" 1 and [32, 35] /3„ ~ 6.0 x 10~ 6 (cm we)" 1 which, 
moreover, reproduces best [5] the Monte Carlo result of Lipari and Stanev [37] for the 
average muon-range in standard rock for E > 10 3 GeV. Furthermore we choose [32] 
the final muon energy to be larger than E mm = 1 GeV. (It should be noted that here 
Pn(E,X) ~ 1 in (27), i.e., R^(E) ~ V max in (28).) The muon-range calculated within 
the continuous muon energy loss approach yields again, as in the case of taus, very similar 
results as the stochastic Monte Carlo calculations [32, 35, 38] for R s *° as shown in Fig. 14. 
This is contrary to the conclusions reached in [10] that the continuous approach to the 
muon energy loss overestimates the muon-range as compared to stochastic Monte Carlo 
simulations. Therefore the continuous approach to the lepton energy loss is applicable to 
both taus and muons, since in both cases it yields similar results for the lepton ranges as 
the stochastic Monte Carlo simulations with the energy loss being treated separately. 
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The transport equation of muon neutrinos includ- 
ing secondary muon neutrinos from tau neutrino 
interactions 



It has been pointed out [14] that the u T — r regeneration chain v T — > r — > v T — > . . . 
creates a secondary + i/^ flux due to the prompt, purely leptonic, tau decays t~ — > 
u T ii~u^ and r + — > z/ r /i + z/^. This will enhance the regenerated fluxes calculated 
according to (1) and thus also the 'naively' calculated [5, 15, 41] upward-going muon 
event rates at the detector site. Secondary neutrinos originate from the associated r 1 * 1 flux 
F T (E,X) and a prompt r-decay like r~ — > v^X'. Adding those contributions, denoted 
by G (_) { } (E,X), to the simple transport equation (1) used thus far one obtains 

with G = G T +^ Ufi where 

G T+ ^ (E, X) = x ^ Ej9) jf * ^ 1/) ^ (E y , X) (30) 

and a similar transport equation holds for F Vji with an appropriate expression for Gr T -_>p . 
The relevant r fluxes F T ± have been calculated in the previous Section (cf. Figs. 10 and 
11). As in (14), the decay kernel in (30) is K^ c (E,y) = (1 - y) dn T +^(z)/dy with 
z = 1 — y and the relevant r + — > v^X' decay distribution is given by [28] 

driT+ ^ z "> = B ^ [ 2 - 6z 2 + Az 3 + P(-2 + 12z - 18z 2 + 8z 3 )] (31) 

with P = +1 and the branching fraction Pv M = 0.18. For a decaying r~ — > ^X' one 
has P = — 1 in (31). Notice that the spectrum in (31) is a little softer than the ^v T 
spectrum from the r 1 * 1 — ^v T X' decay [7, 28] in (15). 

It should be noticed that the contribution of secondary neutrinos may alternatively 
be calculated using directly the V T fluxes F ( _) (E,X) which give rise to the reaction 

V T 
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cc cc 
chains v T — > r~ — > v^X' and v T — > r + — > f M X'. Denoting these contributions by 

G Ur -> D (E, X) and G Dt ->„ (E, X), respectively, the inhomogeneous term in the transport 

equation (29) is given by G = Gp T ^ Ufl with 

where i^/z = E/(l — the decay distribution is given by (31) and the relevant flux 
F Pt has been calculated in the previous Section (cf. Figs. 10 and 11). Although (32) and 
(30) yield the same quantitative results for F v (E,X), these two expressions should not 
be added since it would correspond to double-counting the effect of secondary neutrino 
production. This is due to the fact that the CC contribution G Ut ^ t has been already 
included in (13) [third term on the rhs] for calculating F T . (The situation here is very 
similar to the calculation of the atmospheric muon flux [31, 28] where almost all muons 
come from meson decays with the meson flux being generated by nucleon interactions 
with air, i.e., by nucleon — > meson transitions. These latter transitions are taken into 
account only in the evolution equation of the meson flux, but not anymore for the muon 
flux evolution.) For definiteness, we use the simpler expression in (30) for our subsequent 
calculations. 

As in our previous cases, the most general transport equation (29) for muon neutrinos 
is easily solved by an ansatz like (16) for tau neutrinos, 

X 



F^(E,X) = Fl(E)exp 
with 



Kg(e,x) 



(33) 



and Z v a = Z v + Zq- Inserting (33) into (29) one obtains 

Z Vlt (E,X) = ^j\x< j dyK* c (E,y)r ]l/ (E,y)e- x ' D ^ E > E y^ (35) 
which is similar to (5) but with D Vfi (E, E y , X') = A^ G (E y ,X') - A^ G (E,X'), and 

Z G (E, X) = ±W 1 jf * dX> G(E, X') e x 'l^iEX) . (36) 
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Using again an iteration algorithm to solve for Z v g{E,X), the solution of (35) and (36) 
after the n-th iteration becomes 

rfg\E,X) = l^rfX^ 1 ^^ NC (E,y)^(E, y )e-^< ) ^^) 

^ dX'G(E,X')e X ' /A % G (E,X') (37) 
X Jo 



K{E) 1 



where is as defined above with A^ G ^ A% and A%(E,X') = K{E)/ 

1 — Z^ G (E, X') . Due to the dominant and large r-decay contribution G(E,X) in 
(29), it turns out that the optimal input choice for providing sufficiently convergent iter- 
ative solutions is obtained by implementing, as in the case of tau neutrinos in Sec. 3, the 
peculiar E and X dependence as implied by the r-decays, i.e., by G in (36) from the very 
beginning. Therefore we use again (see discussion after Eq. (20)) zl°J(E,X) = and a 
vanishing Z-factor on the rhs of Z G in (36) which gives for the total input Z-factor 



r(0) 



:(E, X) = ^ 1 jf dX> G(E, X>) e*7M*) . (38) 



Inserting this into the rhs of (37) results in the first iterative solution Z^ G (E,X), and so 
on. In contrast to Z^J in (8), Z^ G does not provide us with a sufficiently accurate final 
result, i.e., the maximum difference between Z^J G and Z^ G is here not always less than 
about 5% for some initial cosmic neutrino fluxes and energies. Therefore we have to carry 
out one further iteration, as in the case of tau neutrinos in Sec. 3, by inserting zjp G into 

(2) 

the rhs of (37) in order to obtain Z K V G (E, X) which turns out to be sufficiently close to 
the final result since \Z® G /Z™ G - 1| < 0.02. 

(1 2) 

Our iterative results for Z y vG are shown in Figs. 15 and 16 together with the appro- 
priate input zf^ G in (38) shown by the dotted curves. In order to illustrate the faster 
iterative convergence for increasing 9 (smaller X), the results for 9 = 50° are presented 
in Fig. 15 as well. The sufficiently accurate results Z^J G (E,X) and the similar expres- 

(2) 

sions for Zy G , when inserted into (33), yield the final total fluxes F u ^ +Pfi (E,X) shown 
in Figs. 17 and 18. The effect and importance of secondary neutrinos is best seen by 
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comparing our results (solid and dashed curves) with the usual ones [3, 5, 7, 15] obtained 
just for primary muon neutrinos (G = in (29)) shown by the dotted curves, which cor- 
respond of course to the results obtained in Sec. 2. Our results in Fig. 17 agree with the 
ones obtained in [9], within the approximations made there. 

(-) ( i } 

The corresponding v M initiated upward-going [i event rate per unit solid angle and 
second is calculated according to 

N^ = N A dE u dyAiE^R^EjT) " , F U ^E V1 X) (39) 

Je™i j ay 

with = (l — y)E l/ and the energy dependent area A{E IJ ) of the considered underground 
detectors is taken as summarized in [5]. The muon-range is given by R^E^, E™ in ) = 
^In ^+ffismln - It describes the range of an energetic muon being produced with energy 
E^ and, as it passes through the Earth loses energy, arrives at the detector with energy 
above £™ n . The energy-loss parameters are taken as at the end of the previous Section, 
i.e., a M = 2 x 1CT 3 GeV (cm we)" 1 and ^ = 6x 1CT 6 GeV (cm we) -1 . The integral over 
the neutrino energy E v was, for defmiteness and better comparison [9], performed up to 
a maximum neutrino energy of 10 8 GeV. The differential ^-dependent fx~ + fi + rates for 
rmin = 1Q 4 GeV and Ef n = 10 5 GeV are shown in Figs. 19 and 20. We also include 
the contributions initiated by the primary + flux, for brevity denoted by — > /i, 
and by the v T + v T flux via v T — > r — > ji and the r _ + r + flux via r^/ias discussed in 
the previous Section. The secondary neutrino contributions to the muon event rates have 
their largest relative contributions obviously at small nadir angles, with an enhancement 
over the primary + v T + r initiated rates of up to 40% for the hard E -1 , AGN-M95 
and TD-SLBY initial fluxes. 

At small nadir angles, however, the event rates are smallest and statistics are low. For 
9 >, 60° the event rates are roughly a factor of (more than) 10 larger and the enhancement 
of the overall z/ M + u T + r initiated muon rates (dashed-dotted curves in Figs. 19 and 20) 
can not be larger than about 15%. These results are more explicitly illustrated in Tables 
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2 and 3 where we present, besides the total nadir-angle-integrated rates, also the ones 
integrated over three typical 6*-intervals. (Remember that this amounts to integrating (39) 
over f„ 27r dip f emax d9 sin 9 = 2n f a max d9 sin 9, with 6L; n = 0° and # max = 90° for the total 

•>V -1 Pmin J fmin 

rates.) Since secondary muon neutrinos contribute significantly to a muon excess only 
at small and medium nadir angles, 9 < 60°, the primary event rates (shown in brackets 
in Tables 2 and 3) can be enhanced by more than 20%, in particular for E™ m = 10 5 
GeV. The statistics, however, are low since the fluxes are already strongly attenuated 
for 9 < 60° (large X), cf. Figs. 17 and 18. On the other hand, most of the events 
are generated at large 9, 9 > 60°, where the effect of secondary neutrinos is sizeably 
reduced (cf. Figs. 19 and 20), the total rates in Tables 2 and 3 are increased by less than 
10%. Since the expected angular resolution of present and proposed detectors [21, 42] is 
typically about 1°/ (i^/TeV) - 7 , differential ^-dependent measurements should be feasible, 
in order to delineate experimentally the effects of secondary neutrino fluxes. Keeping in 
mind that the lifetime of the planned experiments is roughly ten years, it appears to be 
not unreasonable that the tenfold rates implied by Tables 2 and 3 may be observable in 
the not too distant future. 



5 Summary and Conclusions 



For the sake of completeness we have first studied the solutions of the single transport 
equation for cosmic neutrinos propagating through the Earth. Although frequently 
used, the excellent convergence of its iterative solutions has not been explicitly demon- 
strated thus far for more realistic and hard cosmic neutrino fluxes. Using the symbolic 
ansatz for the solution F l/ (E,X) = F®(E)exp [(1 — Z)X/\ V \, with \ u being the neutrino 
interaction length, the most simple input choice zl°J(E,X) = suffices to produce a 
sufficiently accurate iterative result Z^J already after the first iteration, for all presently 
used initial cosmic model fluxes F°. 
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Turning to the iterative solutions of the far more complicated coupled transport equa- 
tions for ^v T and their associated r^ 1 fluxes, a new semi-analytic input algorithm is pre- 
sented which allows for a fast convergence of the iterative solutions: already a second n = 2 
iteration suffices for obtaining a sufficiently accurate result and thus for the final 

V T 

and its associated F T ± fluxes. In order to achieve this one has to implement the peculiar 
E and X dependence as implied by the r ± decay contributions already into the initial 
zeroth order input Z^°\E, X). Choosing a vanishing input Z^ = as in the case for 
fluxes or even the final solution for the ^ flux as an input, Z^> = Z^ , as frequently 
done, results in a far slower convergence of the iterative procedure. For completeness we 
briefly outline also the implications for the upward-going /i~ + fi + event rates for un- 
derground neutrino detectors using some relevant cosmic neutrino fluxes. These events 
are generated by the so called 'primary' \ and fluxes via the weak transitions 

and decays //, v T r — > fi and r — > /i. Furthermore, for calculating the range 

R T (E) of tau-leptons, their energy loss can either be treated 'continuously' by including 
it directly in the transport equation, or 'stochastically' by treating it separately. Both 
approaches give very similar results for R T up to highest energies of 10 12 GeV relevant 
at present. A similar agreement is obtained for the muon range R^(E). Therefore the 
continuous approach is applicable to both taus and muons. This is contrary to claims in 
the literature that the continuous approach overestimates as compared to stochastic 
Monte Carlo simulations. 

Finally, we generalized the single transport equation for by taking into account 
the contributions of secondary and fluxes. These so called 'secondary' muon neu- 
trino fluxes originate from prompt r 1 * 1 decays where the r-leptons are generated by the 
regeneration chain v r — > r — > u T — > . . . when a cosmic v T passes through the Earth. Thus 
the secondary + flux arises from the associated r 1 * 1 flux, as obtained from the coupled 
transport equations for v T and r, which initiates the r — > transitions (t~ — > v T ii~v ii 
and r + — > u T ji + u^). In order to achieve a sufficiently fast convergence of the iterative 
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solutions of the single generalized transport equation of muon neutrinos, one again has 
to implement the peculiar E and X dependence as implied by the weak r-decays already 
into the initial zeroth order input Z^(E, X). In this case one needs only n = 2 iterations 
for obtaining a sufficiently accurate result Z^ 2 \E ) X) for calculating the final secondary 
v M and fluxes. The /i~ + /x + event rates initiated by the secondary neutrinos are largest 
obviously at small nadir angles (9 < 60°), with a relative enhancement of at most 40% 
over the primary z/^ + z/ T + r initiated rates for the hard initial cosmic fluxes like AGN-M95 
and TD-SLBY. At larger nadir angles, 9 £ 60°, the muon rates are dominantly initiated 
by the primary + v T + r flux and the secondary + flux becomes naturally less 
relevant. Thus the secondary neutrino flux will enhance the total nadir-angle-integrated 
muon event rates only by less than 10%. Nevertheless, it should be possible to observe 
the effects of secondary neutrinos with differential ^-dependent measurements, keeping in 
mind that the angular resolutions of the proposed underground neutrino telescopes will 
reach sub-arc-minute precisions. 
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Table 1: Total nadir-angle-integrated upward-going /i" + /i + event rates per year from 
[y T + v T )N and (z/ M + v^N interactions in rock, with the latter being given in parentheses 
which are taken from Table 1 of [5], for various muon energy thresholds E™ m and the appro- 
priate dominant cosmic neutrino fluxes in Fig. 1. The v T + v T and t~ + r + initiated rates are 
calculated according to Eqs. (12) and (14) of Ref. [13], which are added to the + initiated 
rates in parentheses in order to obtain the final total rates. A bar signals that the rates fall 
below 0.01. This table corrects Table I of Ref. [13]. 



Flux 


Detector 


10 3 


Muon-enerj 
10 4 


?y threshold E 
10 5 


min/ GeV 

10 6 


10 7 




ANTARES 


16.63 (13.7) 


6.28 (5.00) 


2.51 (1.98) 


1.06 (0.90) 


0.34 (0.32) 


AGN-M95 


AMANDA-II 


34.90 (29.1) 


10.76 (8.62) 


3.78 (2.98) 


1.58 (1.34) 


0.49 (0.46) 




IceCube 


170.24 (143) 


41.72 (33.7) 


14.22 (11.2) 


5.93 (5.04) 


1.83 (1.74) 




ANTARES 


0.75 (0.60) 


0.39 (0.32) 


0.10 (0.08) 


0.01 (0.01) 




GRB-WB 


AMANDA-II 


1.39 (1.10) 


0.68 (0.56) 


0.15 (0.13) 


0.02 (0.02) 






IceCube 


5.55 (4.35) 


2.59 (2.13) 


0.57 (0.49) 


0.07 (0.06) 






ANTARES 


0.84 (0.62) 


0.59 (0.45) 


0.33 (0.26) 


0.14 (0.12) 


0.05 (0.05) 


TD-SLBY 


AMANDA-II 


1.33 (0.97) 


0.91 (0.68) 


0.49 (0.39) 


0.21 (0.18) 


0.07 (0.07) 




IceCube 


5.11 (3.70) 


3.42 (2.57) 


1.84 (1.47) 


0.78 (0.68) 


0.26 (0.25) 
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Table 2: Nadir-angle-integrated upward-going [i~ + /i + event rates per year for muons with 
energy above E™ m = 10 4 GeV. The events produced by the primary + P M , v T + v T and 
r~ + t + fluxes are given in parentheses, which are obtained from appropriately integrating the 
relevant dashed-dotted curves in Fig. 20. Notice that the total (0° < 9 < 90°) event rates in 
brackets in the last column agree of course with the final total rates in Table 1. Adding to these 
conventional primary rates the ones induced by the secondary + z/ M fluxes, originating from 
t + — > z/ M and r~ — > one obtains the final results shown. These additional secondary + 
contributions are calculated according to (39) and correspond to integrating appropriately the 
relevant solid curves in Fig. 20. 



Flux 


Detector 


0° < 9 < 30° 


Number 
30° < 9 < 60° 


of events 
60° < 9 < 90° 


Total 


AGN-M95 


ANTARES 

AMANDA-II 

IceCube 


0.19 (0.18) 
0.40 (0.38) 
1.62 (1.54) 


1.13 (1.05) 
2.15 (2.03) 
8.54 (8.06) 


5.25 (5.04) 
8.66 (8.34) 
33.35 (32.13) 


6.58 (6.28) 
11.22 (10.76) 
43.50 (41.72) 


GRB-WB 


ANTARES 

AMANDA-II 

IceCube 


0.02 (0.01) 
0.03 (0.03) 
0.12 (0.11) 


0.10 (0.09) 
0.18 (0.17) 
0.69 (0.67) 


0.28 (0.28) 
0.48 (0.48) 
1.84 (1.82) 


0.39 (0.39) 
0.69 (0.68) 
2.65 (2.59) 


TD-SLBY 


ANTARES 

AMANDA-II 

IceCube 


0.01 (0.01) 
0.01 (0.01) 
0.05 (0.04) 


0.07 (0.06) 
0.12 (0.10) 
0.44 (0.38) 


0.54 (0.52) 
0.83 (0.79) 
3.14 (3.00) 


0.62 (0.59) 
0.96 (0.91) 
3.63 (3.42) 
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Table 3: As in Table 2 but for Ef n = 10 5 GeV. A bar signals that the rates fall below 0.01. 



Flux 


Detector 


0° < 9 < 30° 


Number 
30° < 9 < 60° 


of events 
60° < 9 < 90° 


Total 


AGN-M95 


AIM lArtri/o 

AMANDA-II 
IceCube 


U.Ul ^U.Ul ) 

0.02 (0.02) 
0.08 (0.06) 


U.ZU ^U.l / ) 

0.31 (0.26) 
1.16 (0.98) 


3.68 (3.50) 
13.89 (13.20) 


9 (\7 (0 ^1 ^ 
Z.O ( yZ.DL J 

4.01 (3.78) 
15.11 (14.22) 


GRB-WB 


ANTARES 

AMANDA-II 

IceCube 


0.01 (0.01) 


0.02 (0.02) 
0.03 (0.03) 
0.11 (0.10) 


0.08 (0.08) 
0.12 (0.12) 
0.47 (0.46) 


0.10 (0.10) 
0.15 (0.15) 
0.58 (0.57) 


TD-SLBY 


ANTARES 

AMANDA-II 

IceCube 


0.01 (0.01) 


0.03 (0.02) 
0.04 (0.03) 
0.14 (0.12) 


0.32 (0.30) 
0.48 (0.45) 
1.79 (1.71) 


0.34 (0.33) 
0.52 (0.49) 
1.94 (1.84) 



28 



10" 



AGN-SS 




^o' 10° 
E v [GeV] 



10 



12 



Figure 1: Representative differential fluxes of muon neutrinos [y^ + v^) at the production 
site from active galactic nuclei (AGN-SS [22] and AGN-M95 [23]), gamma ray bursts 
(GRW-WB [24]), toplogical defects (TD-SLBY [25] and TD-SLSC [26]) and Z-bursts 
[27]. Due to naive channel counting in pion production and decay at the production site 
[y e : : v T = 1 : 2 : 0) and maximal mixing, u e : : v T — 1 : 1 : 1, these fluxes 
are divided equally between e-, fi- and r-neutrinos when they reach the Earth's surface 
(i.e. will be divided by a factor of 2). The diffuse neutrino flux upper limit of AMANDA 
[20, 21] are shown by the bars with arrows. 
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Figure 2: The ^-factors for neutrinos, as iteratively calculated according to (6), for 
the generic initial neutrino fluxes in (10) and (11) divided by 2. The result for the first 
order iteration is given in (8). For nadir angles 9 > 50°, the second order iterative 
result zj?J becomes almost indistinguishable from Z$ . 
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Figure 3: As in Fig. 2 but for some typical initial cosmic fluxes shown in Fig. 1. The 
(small) TD-SLSC and Z-burst fluxes in Fig. 1 result in a similar iterative convergence 
as the F° ~ E~~ l flux in Fig. 2. 
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Figure 4: The Z-factors for v T neutrinos, Z^ = zj^' + Z\}, as iteratively calculated for 
n — 1, 2, 3 iterations using the input Z^ of (21) which is shown by the dotted curves. 
The generic initial fluxes are taken from (10) and (11) divided by 2. 
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Figure 5: As in Fig. 4 but only for the relevant n — 1, 2 iterations for 9 = 0° and for the 
dominant initial cosmic fluxes in Fig. 1. 
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Figure 6: The ratios of the v T fluxes / F^ according to two consecutive iterations. 

The fluxes are calculated according to (24) using the appropriate Z^> factors shown in 
Fig. 4 for the two generic E^ 1 and E~ 2 initial fluxes. 
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Figure 7: As in Fig. 6 but with the appropriate factors shown in Fig. 5 as obtained 
for the dominant initial cosmic fluxes in Fig. 1. 
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Figure 8: The final total u T + u T fluxes, calculated according to (24) using the appropriate 
iterative n = 2 results Z^ for 9 = 0°, 30° and 60° {Z^ for 9 = 0° is shown in Fig. 4). 
The generic initial fluxes F® t+Dt (E) are given in (10) and (11). 
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Figure 9: As in Fig. 8 but for the dominant initial cosmic fluxes F® t+9t (E) in Fig. 1. For 
9 = 0° the relevant factors are shown in Fig. 5. The TD-SLBY results are multiplied 
by 1CT 1 as indicated. 
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Figure 10: Attenuated and regenerated v T + v T and r~ + r + fluxes calculated according 
to (24) and (25), respectively, using the sufficiently accurate second iterative solution 
Z( 2 \E, X) for nadir angles 9 = 0° and 30°, and the generic initial F® t+Pt (E) flux in (10) 
which is shown by the dotted curve. The results for the + r + fluxes are multiplied by 
10 5 as indicated. 
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Figure 11: As in Fig. 10 but for the generic E 2 flux in (11) and the relevant dominant 
initial fluxes F® t+Pt = \ d^/dE v in Fig. 1 which are shown by the dotted curves. 
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Figure 12: The r-decay length A^ ec = (E/m T )cT T p, shown by the dashed curve, and the 
r-ranges in standard rock (p = 2.65 g/cm 3 ) for an incident r-energy E and final r-energy 
larger than E mm = 50 GeV. R T is calculated according to (27) and (28); omitting A^ c 
in (27) results in i? T (no A^ c ). The stochastic Monte Carlo evaluations of Rf° are shown 
by the upper dotted curve [38] and the lower one according to [32, 35] which are based 
on the ALLM97 parametrization of structure functions for calculating the photonuclear 
cross section. Using a different parametrization (Bugaev-Schlepin) for extrapolating the 
latter cross section to ultrahigh energies, results in a r-range which is about 25% smaller 
than the upper dotted curve at E — 10 12 GeV [38]. 
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Figure 13: The CC r-interaction length A^? c (upper solid curve) and the r-decay length 
neglecting the energy-loss (dashed curve) and including the energy-loss in standard rock 
(lower solid line). The latter two curves correspond to the ones shown in Fig. 12. 
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Figure 14: The /i-ranges in standard rock (p = 2.65 g/cm 3 ) for an incident muon energy 
E and the final //-energy larger than E mm = 1 GeV. is calculated according to (27) 
and (28) as discussed in the text. The stochastic Monte Carlo evaluations of R s ^° are 
taken from [32, 35] (lower dotted curve) and from [38] (upper dotted curve). 
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Figure 15: The Z-factors for primary and secondary neutrinos, Z V J G = Z V J + Z G , as 
defined in (34), (35), (36) and calculated iteratively for n — 1, 2 according to (37) using 
the input Z^ G of (38) which is shown by the dotted curves. The generic initial E^ 1 flux 
is taken from (10) divided by 2. 
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Figure 16: As in Fig. 15 but only for 9 = 0° and for the generic initial E 2 flux 
divided by 2, and the three dominant initial cosmic fluxes in Fig. 1. 
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Figure 17: The final total + fluxes for primary and secondary muon neutrinos, 
calculated according to (33) using the appropriate iterative n = 2 results for Z^ G for 
8 = 0° and 60°. The generic initial fluxes F° +P (E) are given in (10) and (11). The 
usual primary + fluxes are for brevity denoted by (dotted curves). Similarly the 
secondary neutrino contributions due r + — > and t~ — > are denoted by r — > z/ M . 
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Figure 19: The + fi + event rates for muons with energy above 10 4 and 10 5 GeV 
originating from (anti) neutrinos with E v < 10 8 GeV for the generic initial E^ 1 flux in 
(10). The muon events initiated by the primary and {y^+u fl ) + {u T +u T ) + {T^ +t + ) 

fluxes are for simplicity denoted by and + u T + r, respectively (dotted and dashed- 
dotted curves). The additional events initiated by the secondary muon neutrinos arising 
from r~ — > r + — > z/ M are for brevity denoted by (r — > z/ M ). 
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